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Ion-implantation is a useful technique to study irradiation damage in nuclear materials. To study He effects in 
nuclear fusion conditions, He is co-implanted with damage ions to reproduce the correct He/dpa ratios in the desired 
or available depth range. However, the short-term fate of these He ions, i.e. over the time scales of their own collisional 
phase, has not been yet unequivocally established. Here we present an atomistic study of the short-term evolution 
of He implantation in an Fe substrate to approximate the conditions encountered in dual ion-implantation studies in 
ferritic materials. Specifically, we calculate the fraction of He atoms that end up in substitutional sites shortly after 
implantation, i.e. before they contribute to long-term miscrostructural evolution. We find that fractions of at most 3% 
1 should be expected for most implantation studies. In addition, we carry out an exhaustive calculation of interstitial 
He migration energy barriers in the vicinity of matrix vacancies and find that they vary from approximately 20 to 60 
meV depending on the separation and orientation of the He-vacancy pair. 
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This paper tries to answer a simple yet important ques- 
tion in He-implantation studies: Do ion-implanted He 
atoms end up as interstitial or as substitutional particles 
in the target matrix? The difference is critical because 
of the large diffusivity difference between both forms of 
He: interstitial He (i-He) diffuses extremely fast, sam- 
pling large portions of the configurational space quickly, 
readily finding other defects or microstructural features. 
Conversely, substitutional He (s-He), while energetically 
more stable, is immobile, necessitating migration of other 
point defects before it can move. This can happen either 
by reacting with a self-interstitial atom (SIA) that recom- 
bines with the vacancy and knocks the He back to an in- 
terstitial site, or by correlated lattice exchange reactions 
with a vacancy in nearest-neighbor positions. Either way, 
the diffusivity of s-He is still several orders of magnitude 
lower than that of i-He. 

\ Despite the important implications of these mecha- 
nisms on the subsequent microstructural evolution, at 
present most researchers consider that 100% of the im- 
planted He is interstitial and can only become 
substitutional by finding an isolated vacancy in an un- 
corrected fashion via long-range diffusion. The question 
we ask here is whether this is true for all He atoms or 
whether some of them can become substitutional as part 
of their own implantation process prior to uncorrelated 
diffusion taking place. 

Here we present a computational study involving the 



binary collision approximation (BCA), molecular dynam- 
ics (MD), and kinetic Monte Carlo (kMC) simulations. 
The BCA is used to simulate the penetration of He beams 
of various energies into Fe targets, and to obtain en- 
ergy distributions of Fe recoils due to He impact. MD 
is then used to simulate He thermalization after the pri- 
mary knock-on event in its own collisional environment, 
and to ascertain whether He atoms create stable Frenkel 
pairs that can result in correlated recombination. Finally, 
we use kMC to calculate the fraction of freely-migrating 
He atoms from those that do create defects but do not 
find the vacancy during MD time scales. From these sim- 
ulations, we find that nearly 3% of the He atoms end up 
in subsitutional sites. While this number appears small, 
it nonetheless leads to dramatic differences in the mi- 
crostructural evolution of the material, as will be shown 
in future studies. 

2. Results 

2.1. Calculation of recoil distributions and He energies 

The He energy range of interest for fusion materials 
lies between 3.5 and 0.33 MeV, corresponding to the en- 
ergy of a particles emitted from fusion reactions and 
those produced via (n,a) transmutation reactions in Fe 1 . 
In addition, ion beam experiments typically use energies 
in this range to achieve penetrations of a few microns, 
so it is useful to have an intermediate energy for refer- 
ence. Thus, we first calculate the Fe recoil distribution 
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Figure 1: Cumulative Fe recoil distribution for He-ion beam irradia- 
tions with recoil energies of 0.33, 1.7 and 3.5 MeV incident energy. The 
threshold displacement energy is 25 eV. 



Table 1: SRIM parameters for the three He-beam energies considered. 



He ion energy (MeV) 


0.33 


1.7 


3.5 


Depth (fim) 


0.7 


2.6 


6.0 


% energy to recoils 


0.43 


0.11 


0.07 


Average recoil energy (eV) 


194 


211 


222 


Maximum recoil energy (keV) 


66 


134 


315 



for three He-ion energies, namely 0.33, 1.7 and 3.5 MeV, 
using SRIM yfl. The cumulative recoil energy distribu- 
tions in each case are given in Fig. [TJ where a threshold 
displacement energy of 25 eV was used. As the figure 
shows, the three recoil distributions are almost identical. 
This is because He ions only create recoils when they 
have slowed down to a few keV, without much participa- 
tion from their higher-energy histories, which as shown 
in Tab.[T]only contribute to penetration and maximum re- 
coil energy. It is more informative to compare the average 
recoil energies, which, in contrast, differ only by a few eV. 
In all cases, the energy expended in recoils amounts to 
less than 0.5% of the total ion energy. 

These results suggest that in the energy range relevant 
to fusion materials the actual ion energy is irrelevant for 
damage purposes. Therefore, we take the 1.7 MeV spec- 
trum shown in Fig. [2] as representative of all He energies 
and proceed to simulate the effect of these recoils on lat- 
tice damage. We note that SRIM does not capture chan- 
neling, which may have some impact on the final results. 

2.2. Molecular dynamics simulations of He impact in Fe 

Next we study the fate of He ions at the end of their 
collision trajectories, when they collide with the last of 
their recoils before thermalizing in the host lattice. We 
assume that these recoils are ejected with energies consis- 
tent with the recoil spectrum obtained using SRIM (see 
Fig. |2J, as the probability for high energy recoils, pro- 
duced early in the He collision sequences, is very small. 



Figure 2: Spectrum of Fe recoils sampled by MD simulations in com- 
parison with spectrum for 1.7 MeV He obtained using SRIM. 



The objective is then to investigate whether He ions can 
become substitutional by interacting with vacancies of 
their own creation, rather than by long-range diffusion. 
For this, we perform MD simulations of He collisions in 
Fe and analyze the final configurations. 

Simulations were carried out using the massively par- 
allel MD code lammfs The simulation cell consisted 
of 38 x 38 x 38 conventional body-centered cubic (BCC) 
unit cells corresponding to 109,744 Fe atoms. All simula- 
tions were carried out at a temperature of 700 K, which is 
representative of fusion conditions. Atomic interactions 
were modeled using the Fe-He interatomic potential of 
developed by Juslin and Nordlund J/J, which builds on 
the Fe potential by Mendelev et al. [8j and gives an equi- 
librium lattice constant of a = 2.871 A at 700 K. The col- 
lision simulations were done with He as an energetic ion 
in an otherwise perfect Fe lattice. The He atoms was 
assigned kinetic energies that produce a Fe recoil dis- 
tribution consistent with the SRIM data in Fig. [2] We 
sample the Fe recoil velocity upe directly from the SRIM 
spectrum and, assuming purely elastic collisions, assign 
a velocity vyi e to the He atom 
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where the mass ratio is m^i e / mp e « 0.077 and the velocity 
vector is directed at the Fe PKA. 

The atoms that were initially within a radius of 18.65 a$ 
of the center of the simulation cell were considered the 
core region while all other atoms were assigned to the 
edge region. The latter were coupled to a Nose-Hoover 
thermostat at 700 K, while the remainder of the system 
evolved in time according to the microcanonical ensam- 
ble. The equations of motion were integrated until ei- 
ther the He atom thermalized (kinetic energy less than 
0.1 eV) or escaped the core region. The final configura- 
tion was relaxed using conjugate gradient minimization 
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Figure 3: Average number of stable Frenkel pairs created during the 
MD simulations of He collision cascades in Fe. The bars indicate the 
standard deviation of the number of defects created. 

and defects were identified by analyzing the occupancies 
of the Wigner-Seitz cells of the initial perfect BCC lattice. 
In total, we simulated 5,000 events to extract sufficient 
statistics. 

Out of the total 5,000 events simulated, 2,668 (53.4%) 
were terminated when the He atom exited the core re- 
gion. They are interpreted as cases in which the He equi- 
librates in a region "far" away from the PKA, such that 
the probability for the He to recombine with defects cre- 
ated in the original cascade is very small. Accordingly 
these events are counted as i-He in the total atom tally. 
The remaining 2,332 events were analyzed to obtain the 
number and distribution of irradation induced point de- 
fects. In 25 cases the helium ended up in a substitutional 
site corresponding to 1.07% of the "thermalized" cases 
and 0.5% of the total number of events. This occurred 
typically within 4 to 6ps of simulated time. Of the re- 
maining events, 673 cases (13.5%) resulted in He ther- 
malization but no Frenkel pairs, further adding to the 
i-He tally 

The remaining 1,634 cases (32.7%) deserve special at- 
tention. This group includes those He atoms that have 
thermalized within the core region and have created sta- 
ble Frenkel pairs but have not become substitutional dur- 
ing the MD simulation. The number of stable point de- 
fects created by the He collisions as a function of PKA 
energy is shown in Fig. [3] Although in principle, He can 
also be trapped by SIAs J§], we did not observe any in- 
stances where this occurred. SIAs created during the cas- 
cade either diffused away from the core region or were 
situated too far from the final position of the He atom. 

Analysis of the resulting configurations enables us to 
determine the spatial correlation of He interstitials with 
vacancies as shown in Fig. [4] The overall distribution has 
a mean of approximately 30 A, with lower energy con- 
tributions slightly shifted to shorter separations. In any 
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Figure 4: Pair correlation of He interstitials and Fe vacancies obtained 
from 1,634 cases of He impact in a Fe lattice. Comparison of the relative 
contributions of different energy windows shows a slight shift toward 
smaller distances for smaller energies. 

case, to determine the fate of these He ions conclusively, 
one needs to "age" these configurations further using a 
technique capable of probing longer time scales. This is 
akin to calculating the fraction of freely migrating defects 
in high-energy cascade simulations (lOfl . To this end, we 
carry out kinetic Monte Carlo (kMC) simulations of He- 
vacancy reactions according to the distribution given in 
Fig. II 

2.3. Kinetic Monte Carlo simulations of He-vacancy reactions 

The kMC simulations consisted of a single vacancy lo- 
cated at the center of a reaction sphere and a randomly- 
oriented He atom separated from the vacancy by a dis- 
tance sampled from the distribution shown in Fig. |4] We 
have used the Green's function Monte Carlo method ifllll 
in a continuum Fe medium with two spherical particles 
representing the He atom and the vacancy. The sum of 
the radii of the vacancy and He atom was set equal to 
the third-nearest neighbor distance in the BCC lattice, 
r = an a/2 to be consistent with the He-V binding en- 
ergy calculations performed in the Appendix. Other au- 
thors have suggested that binding occurs up to the fifth 
nearest neighbor distance ifTill . No further correlation 
between the He atom and the vacancy is assumed, i.e. 
jumps toward and away from the central vacancy are 
sampled with equal probability. The critical parameters 
for the kMC simulations are the temperature, the diffu- 
sivities, and the size of the simulation box. The tempera- 
ture was 700 K, the same as in the MD simulations. With 
respect to the diffusivities, we neglect vacancy diffusion, 
as its diffusion coefficient is known to be several orders 
of magnitude lower than that of He, even at these tem- 
peratures. The diffusion coefficient of He in BCC Fe was 
taken from Terentyev et al. who used the same in- 

teratomic potential as in the present work and obtained 
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Table 2: Summary of results from MD and kMC simulations of He impacts in BCC Fe. For the data obtained from kMC simulations the third 
column states the number of occurences equivalent to the number of events treated by MD simulations. The actual number of events simulated by 
kMC was 10 6 (see text for details). 

Event Number of occurences Relative fraction 



i-He He escaped core region during cascade (MD) 2668 53.3% 

no stable defects created during cascade (MD) 673 13.5% 

He escaped core region after thermalization (kMC) 1513 30.3% 

total 97.1% 

s-He He trapped during cascade (MD) 25 0.5% 

He trapped after thermalization (kMC) (MD) 121 2.4% 

total 2.9% 



Dpie — 5.1 x 10~ 3 exp ( — 76/kT) cm 2 s _1 (migration en- 
ergy in meV). Finally, also for consistency with the MD 
simulations, we considered a spherical region with ra- 
dius equal to 18.65flo- 

If over the course of a simulation the He atom escaped 
the spherical region, it was added to the i-He count, 
whereas, if it reacted with the vacancy in the center of 
the simulation sphere, it was tallied as a s-He. After 
10 6 events, we calculated the probability for He-vacancy 
recombination under these conditions to be 7.4%. The 
probability of He escape in an infinite medium has a 
known analytical solution given by 11311 

P = l-| (2) 

where d is the initial separation distance. This formula 
yields an average reaction probability of approximately 
12%, which is higher than the kMC value because it is 
not limited to a finite reaction volume. 

When prorated to the number of cases that consti- 
tute the He-vacancy pair correlation in Fig. [4) this re- 
sult, added to the 0.5% computed directly from the MD 
simulations, results in a total fraction of s-He of 2.9% in 
ion implanted BCC Fe. The various contributions to this 
number are summarized in Tab. [2] for clarity. 

3. Discussion and conclusions 

Understanding how implanted or transmutation He 
is partitioned between i-He and s-He is paramount be- 
cause of the different characteristics of both species in the 
BCC lattice. A comprehensive review on the role of He 
in Fe has been recently published [14], where the basic 
energetics are given and a review of the existing litera- 
ture is provided. All studies agree that s-He is energeti- 
cally more favorable with a lower formation energy and a 
strong He-vacancy binding energy but cannot move un- 
less aided by other point defects. For its part, i-He dif- 
fuses very fast in all three dimensions, rapidly probing 
vast regions of configuration space and finding sinks or 



other defects very efficiently. These different behaviors 
have important implications in terms of the long term 
microstructural evolution. For example, He in solution in 
the BCC lattice is known to stabilize vacancy clusters pro- 
duced directly in high-energy cascades. In terms of its ef- 
fect on direct damage production, however, the evidence 
reported in the literature is contradicting. On the one 
hand, some researchers using the Fe-He Wilson-Johnson 
potential [15] have found that high-energy cascades in 
BCC Fe doped with small concentrations of s-He result 
in higher numbers of vacancy clusters than in pure Fe 
[16]. They also found these clusters to be generally larger 
in size. On the other hand, using the Juslin-Nordlund 
potential — employed here — Lucas and Schaublin have 
found that it is i-He in solution that causes larger clus- 
ter sizes and number densities to appear [17]. In fact, 
they observe that s-He reduces the number of stable de- 
fects with respect to pure Fe. These workers also per- 
formed a systematic Fe-He potential comparison, not- 
ing that the Juslin potential gives results that are overall 
in better agreement with DFT calculations 1 18]. In any 
case, these results show the importance of determining 
the correct partition of implanted He, something typi- 
cally neglected in most rate theory studies, where He is 
generally inserted as i-He (although some notable excep- 
tions exist 

B). 

Next we discuss the validity and limitations of our 
approach. This work hinges on the fact that Fe recoil 
spectra from He ions with a wide range of incident ener- 
gies are almost identical. Further, assuming that recoils 
produce isolated cascades, data extracted from ion beam 
experiments can be used to infer the behavior of a parti- 
cles created in (n,oc) reactions. Then, the only difference 
between (n,a) reactions and He-ion irradiations is that 
the former are created homogeneously within the matrix, 
whereas ions penetrate a short distance into the material 
(c.f. Table [I), but up to the cooling-down phase of the 
cascade, damage is produced in an identical fashion in 
both instances. 

With regard to the interatomic potentials used, we 
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have already mentioned the studies in Refs. ifl^. 
Stewart et al. jioll have recently carried out an exhaustive 
comparison of Fe-He and He-He potentials available in 
the literature. These authors noted that the Fe-He poten- 
tial is strongly dependent on the matrix (Fe-Fe) to which 
it is coupled, and that the potentials used in the present 
work produce little clustering. In addition, Yang et al. 
lElll and Pu et al. fzj have shown that different poten- 
tials can have a noticeable influence on vacancy cluster 
formation in displacement cascades. All of these effects 
should again be mitigated at high temperatures, yet in 
light of these results some variability in the final fraction 
of s-He can be expected. 

Another, limitation is the assumption of uncorrelated 
diffusion in the kMC calculations. Presumably, the mi- 
gration energy of an i-He varies as a function of its 
proximity to a vacancy from the value in the bulk (here 
E m = 58 meV, see Fig. |A.5t to zero in the two intersti- 
tial sites closest to a vacancy (leading to spontaneous re- 
combination). In this work, however, we have neglected 
this dependency, for two main reasons. First, because 
at the simulation temperature, the i-He diffusion barri- 
ers, which are typically on the order of 40 to 60 meV (see 
Tab. IA.3t , are well below kinetic energy fluctuations, as 
shown in the analysis provided in the Appendix. Sec- 
ond, this correlated effect is partially captured by adjust- 
ing the interaction distance, which was set to 5 th nearest 
neighbor distance in this work. Indeed, other (lattice) 
kMC calculations of He and He-vacancy complex diffu- 
sion have also disregarded this effect for similar reasons 

In conclusion, we have performed a computational 
study of He implantation in BCC Fe and have calculated 
the fraction of implanted He that becomes substitutional 
during its own collisional phase. We find that a fraction 
of approximately 3% is reasonable for the energy range 
of interest to fusion materials. Damage accumulation cal- 
culations should take this datum into account, although 
more calculations are needed to accuratelt quantify the 
impact on long-term microstructural evolution. 
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Figure A.5: Conventional unit cell of the body-centered cubic lattice 
illustrating the locations of tetrahedral and octahedral interstitial sites, 
and He migration barriers in the defect-free iron. Note that for clarity 
only a subset of the interstitial sites is shown. 

Notes 

1 The transmutation reaction 56 Fe(n,a) 53 Cr results in an excess 
mass of m e = m Fe + m n - (m Cr + m«) = 55.9349375 + 1.0086692 - 
(52.9406494 + 4.0026032) = 0.0003541 amu. This is equivalent to E = 
m e c 2 = 0.33 MeV, although some variability in the form of a relatively 
narrow energy spectrum is to be expected depending on local condi- 
tions such as orientation, atomic vicinity, etc. (Source of particle rest 
masses: NIST fl) 

Appendix A. Helium interstitial migration in the 
vicinity of an iron vacancy 

To supplement the kMC simulations conducated as 
part of the present work, we carried out a systematic 
study of He interstitialimigration in the vicinity of a va- 
cancy. Calculations of migration barriers for He intersti- 
tials in defect-free iron and in the vicinity of a vacancy 
were carried out using the drag method implemented by 
the authors in the MD code lammps. All barrier calcu- 
lations were carried out at the zero-K lattice constant of 
2.855 A using 16 x 16 x 16 supercells based on the con- 
ventional BCC unit cell. 

The formation energies of interstitial helium in tetrahe- 
dral and octahedral sites are 4.39 eV and 4.52 eV, respec- 
tively. In the ideal lattice the coordinates of the tetrahe- 
dral and octahedral interstitial are and (0,0,^) 
in units of the lattice constant. As can be deduced from 
Fig. IA.51 there are 12 tetrahedral and 6 octahedral sites in 
the conventional BCC unit cell. Figure IA.5I also summa- 
rizes the results for the migration barriers of tetrahedral 
He interstitials in defect-free iron. We obtain a barrier 
of 58meV for migration along (110) in agreement with 
Ref. [7]. For migration along (100) the calculations yield 
a value of 127 meV, which precisely corresponds to the 
enegry difference between the tetrahedral and octahedral 
sites as the latter coincides with the saddle point. 
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Table A.3: Migration barriers for He interstitials in the vicinity of a 
vacancy. x m (. initial position in fractional coordinates with respect to 
the vacancy located at the origin, N„„: initial neighbor shell, Xfj„: final 
position, Nji„: final neighbor shell, AE;_f = Ej — £,: energy difference 
(meV), AE;,: migration barrier (meV). 
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Before we address the migration of He in the vicinity 
of a vacancy, we first study the effect of the vacancy strain 
field on the energetics of He interstitials. To this end, we 
constructed all crystallographically distinct He-vacancy 
pairs that can occur within a radius of four lattice con- 
stants, which yields 78 distinct pairs. 

The binding energy for a He-vacancy pair is defined as 

AE b = AE / (He - V) - AE / (He) - AE / (V), (A.l) 

where AE / (He - V), AE f (He) = 4.391 eV, and AE / (V) = 
1.721 eV are the formation energies of the He-vacancy 
pair, the isolated He interstitial as well as the isolated 
vacancy. Negative and positive values of AE;, indicate 
attraction and repulsion of the He-vacancy pair, respec- 
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Figure A.6: Binding energy of He-vacancy pairs as a function of sepa- 
ration after relaxation. 

tively. The formation energy is given by 

AE f = E def - ^E id , (A.2) 

where E, and N, are the total energy and the number of 
Fe atoms in configuration i. Here, we have quietly set 
the chemical potential of He to zero which has a mini- 
mal effect on the formation energies and no effect on the 
binding energy. 

The binding energy is shown as a function of the He- 
vacancy separation in Fig. IA.6I Note that the two nearest 
He-vacancy pairs have been omitted in Fig. IA.61 since 
they spontaneously recombine with the vacancy leading 
to a strongly negative binding energy of —2.01 eV. We 
find notable variations in the binding energy at least up 
to a separation of about 2.2 lattice constants. It is fur- 
thermore noteworthy that the variation in the binding 
enegry is not monotonic with the distance but exhibits 
pronounced oscillations that to some extent can be corre- 
lated with different crystallographic directions. 

We can now consider the different possibilities for He 
interstitials migrating in the vicinity of a vacancy. From 
Fig. IA.51 we can deduce that from any tetrahedral site 
there are six possible jumps, four of which are along 
(110) and two of which are along (100). More specifi- 
cally, 

1. tet -> tet: (0, \, \) + {0,+\, +\) -> (0, \, f ) 

2. tet -> tet: (0, \,\) + (0,+|, -\) -> 

3. tet -> tet: (0, \, \) + {+\, -|,0) ->■ 

4. tet -> tet: (0, \, \) + {-\, -|,0) -> 

5. viaoct: (0, \, \) + (0, +^,0) -3- (0,f,±) 

6. viaoct: (0, \, \) + (0, -^0) ->■ 

The barriers that are obtained for these paths start- 
ing from different initial He-vacancy arrangements are 
summarized in Table IA.3I We find that while there is a 
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great variability of the values, in general migration bar- 
riers range from about 20 to 50meV for shells 3 and 4, 
and from 40 to 60 meV for shells 5 to 7, relatively quickly 
closing in on the bulk value of 58 meV. Using these data, 
we rationalized (see Sec. [3) that a good choice for the 
interaction radius entering the kMC simulations is five 
lattice constants. 
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